Geometry of frictionless and frictional sphere packings 
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We study static packings of frictionless and frictional spheres in three dimensions, obtained via 
molecular dynamics simulations, in which we vary particle hardness, friction coefficient, and coef- 
ficient of restitution. Although frictionless packings of hard-spheres are always isostatic (with six 
contacts) regardless of construction history and restitution coefficient, frictional packings achieve 
a multitude of hyperstatic packings that depend on system parameters and construction history. 
Instead of immediately dropping to four, the coordination number reduces smoothly from 2: = 6 as 
the friction coefficient jj, between two particles is increased. 
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I. INTRODUCTION 

Dense amorphous packings of frictionless spheres have 
proven to be an extremely useful paradigm in different 
physical contexts, such as metallic glasses colloidal 
crystals |^ , and emulsion rheology ^ . Granular materi- 
als are another example of a system with macroscopically 
large particles, with one major difference: grain-grain in- 
teractions involve frictional forces. As a result, granular 
packings may be quite different from frictionless sphere 
packings in ways which may impact significantly on their 
physical properties. 

A common quantity of interest in packings of hard 
spheres is the average number of contacts per particle 
(coordination number) z. In order to achieve static me- 
chanical equilibrium, each sphere in the packing needs a 
sufficient number of constraints that freeze out its trans- 
lational and rotational degrees of freedom. These con- 
straints are provided by contacts, and once there are a 
sufficient number of them, the packing can accommo- 
date external body or boundary forces, as long as a set 
of contact forces satisfying mechanical equilibrium can be 
found for the given arrangement of such contacts. The 
minimal average coordination number required to obtain 
static packings of c?-dimensional frictionless spheres that 
are stable against external perturbations is 2:„ = 2(i [Q, 
whereas for spheres with friction, Zf = d+l In three 
dimensions, z„ = 6 and Zf = A. We call such packings 
"isostatic" . 

In this study, we investigate whether or not sphere 
packings readily achieve isostaticity under generic pack- 
ing conditions. This isostaticity hypothesis is impor- 
tant in theories focusing on the macroscopic response of 
such packings ||^,and at first appears reasonable, given 
the strong numerical evidence from simulation that 
Zn = 6 for dense random packings of frictionless spheres. 
Simulation studies of frictional spheres compressed in a 
gravity- free environment Q have shown that Zf is signif- 
icantly less than 6, but with the lowest achieved value 
of around 4.5, it remains unclear whether the minimal 
value of 4 is reached in the limit of zero confining pressure 



(equivalent to the hard-sphere limit.) It is also unclear 
in what way the packings would change for arbitrarily 
small friction coefficient fi between the spheres in order 
to achieve a reduction in z to four, if the isostaticity hy- 
pothesis were true. 

To address such questions, we perform a systematic 
simulation study of the effect of various parameters on 
sphere packings. In particular, we vary the following ma- 
terials properties: the sphere hardness fc„, the coefficient 
of restitution e, and the friction coefficient /i between 
two particles. We also vary the initial conditions of the 
packing by varying the initial packing density (f)^ , as well 
as the initial velocities of the spheres. We investigate 
how the density, coordination number and the nature of 
the contacts change as these parameters are varied. Al- 
though frictionless hard-spheres appear to form isostatic 
packings regardless of construction history and restitu- 
tion coefficient, frictional hard-spheres achieve a multi- 
tude of hyperstatic packings {z > Zf) that depend on 
system parameters and construction history The co- 
ordination number reduces smoothly from 2: = 6 as the 
friction coefficient is increased, disagreeing with the iso- 
staticity hypothesis. 



II. SIMULATION METHOD 

We present molecular dynamics (MD) simulations in 
three dimensions on model systems of TV = 20, 000 
monodisperse, cohesionless spheres of diameter d and 
mass 771. The system is spatially periodic in the 
xy— plane, with a unit cell of size 20d x 20d, and is 
bounded in the z— direction by a rough bed at the bottom 
and an open top. The starting configurations consist of 
randomly positioned non-overlapping spheres, with pack- 
ing fractions in the range 0.02 < 0' < 0.3, obtained by 
varying the overall height while keeping N and the x, y 
dimensions fixed. The system is subsequently allowed to 
settle under gravity on top of the rough bed (see Fig. g). 
The equilibrated static packing height is about 50d. This 
method of construction mimics the pouring of granular 
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materials through a sieve to an area far away from side 
waUs, without forming a conical heap. 

Different ways of preparing static granular packings 
include compressing them Q, and reducing the inclina- 
tion angle of gravity-driven chute flows below the angle 
of repose In the frictionless case, conjugate gradient 
methods have been used to study dense random packings 
However, for the case of particles with friction, MD 
simulations are more appropriate in order to properly ac- 
count for both the normal and tangential forces, since the 
latter depend on the loading history of the contact. 




(a) (b) (c) 

FIG. 1. Lower portion of the packing of TV = 20,000 
spheres in a periodic cell 20d x 20d supported by a rough 
bed (black particles), constructed by settling under gravity, 
with fi = 0.50, fc„ = 2 ■ W^mg/d, and e = 0.88. (a) Initial 
configuration with volume fraction ~ 0.13.(b) Intermedi- 
ate time during settling, (c) Final (static) configuration with 
(j)-^ 0.60. The black frame is added as a guide to the eye. 

The spheres interact only on contact through a lin- 
ear spring-dashpot interaction law |ic| ] in the normal 
and tangential directions to their lines of centers ||ll|,|l2| . 
Contacting spheres i and j positioned at and expe- 
rience a relative normal compression 5 ~ \rij — d\, where 
Tij = Ti — Yj , which results in a force 

Fy=F„ + Ft. (I) 

The normal and tangential contact forces are given by 

Tfl 

F„ = knSnij - — 7„v„, (2) 



Ft = -fc^Ast - -7(Vt, (3) 



where = Vij/rij, with = |ry|, v„ and Vt are the 
normal and tangential components of the relative surface 
velocity, and kn^t and ^n,t are elastic and viscoelastic con- 
stants respectively. Asj is the elastic tangential displace- 
ment between spheres, obtained by integrating surface 
relative velocities during elastic deformation of the con- 
tact. The magnitude of Ast is truncated as necessary to 
satisfy a local Coulomb yield criterion Ft < /ii^n, where 
Ft = |Ft| and Fn = |F„|. Frictionless spheres can be 
simulated simply by setting ^ = 0. Finally, the parti- 
cles are moved according to the total forces and torques 
applied to them through contacts and the gravitational 
field. Additional details can be found in Ref. p3[ . 

The presented simulations were carried out for a range 
of materials parameters. The normal spring constant kn 
varied from 2 ■ 10^ to 2 ■ IQ^mg/d in order to understand 
how the system behaves as it approaches the hard sphere 
limit. In all cases, kt — 2kn/l |0. In order to study 
the crossover from frictionless to frictional systems, the 
local particle friction coefficient /it was varied from to 
10. Finally the coefficient of restitution e, i.e. the ratio of 
the final to initial normal velocities in a head-on binary 
collision, was also varied. For a linear spring-dashpot 
interaction, 

e = exp(-7„icoi/2), (4) 

where the collision time icoi, 

t,oX=^{2kn/m--il/A)-^'^. (5) 

Three values of e=0.26, 0.50, and 0.88 were used. The 
effect of e is like that of varying the quench rate in a ther- 
mal system: for smaller values of e, collisions are more 
inelastic, energy is dissipated faster, and the spheres have 
less ability to move from the the first position where they 
are mechanically stable, which may result in packings 
that are less stable than similarly prepared packings of 
more elastic grains. 

Most of the simulations were started from a static 
configuration of 0* w 0.13 as shown in Fig. |l|(a). The 
timestep 5t ~ 0.05-\/m//c„ was chosen to accommodate 
the decreasing collision time as the particle hardness 
kn is increased [cf. Eq. (H)]. For fc„ = 2 • lO^mg/d, 
St « lQ~'^yJd/ g. Simulations were then run until the ki- 
netic energy per particle was less than 10~^mgd for small 
kn, and up to three orders of magnitude less for large fc„. 
This requires 3 - 8 • lO'^St for small fc„, and 4 - 8 • W6t 
for fc„ = 2 • lO^mg/d. 

Because of the large amount of time required to reach 
static equilibrium, each set of parameters was run for 1 
to 5 configurations. In some cases, the particles were 
given an initial, random velocity with an average KE per 
particle of approximately 20 — lOOmgd. The results were 
identical in all cases within sample to sample fluctua- 
tions. 
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III. RESULTS 



A. Coordination Number 



Figure || shows the effect of friction coefficient on 2; 
for fc„ = 2 • lO^mg/d and e = 0.88 and 0.26. For both 
values oie, z — 6.1 44±0 .002 for frictionless packings. As 
we shall see in Sec. [II C| , the deviation from the isostatic 
value of 6 can be attributed to the finite stiffness of the 
spheres, the isostatic value is apparently obtained in the 
hard-sphere limit. However, there is no sudden drop from 
z = 6 as friction is turned on; rather there is a gradual 
decrease in z to a parameter dependent minimal value, 
accompanied by a similar decrease in the final volume 
fraction (j)^ (see Fig. |[) 




FIG. 2. Bulk averaged coordination number 2 as a function 
of for fc„ = 2 ■ Wmg/d and (p^ — 0.13 for two values of e. 



As depicted in Fig. ^ the decrease in z is primarily 
due to an overall shift in the distribution of coordina- 
tion numbers to lower values, rather than a change in its 
shape and width. Consequently, the frequency of parti- 
cles with eight or more neighbours reduces as fi increases, 
and particles with as few as two contacts start to appear 
at /i = 0.5, indicative of arching within the packing at 
large fi. The saturation of 2 and (j>-^ for fi ^ 1 is due to 
the fact that the typical tangential forces Ft in a packing 
with /i = cxD is expected to be of order Fn, and lower- 
ing /i has little effect on the packing down to /i « 1. 
This behavior of the contact forces is further verified in 



Sec. HID 



Unlike the frictionless case, the deviations from iso- 
staticity {z = 4) for /i > cannot be attributed to cor- 
rections due to the finite stiffness of the spheres. The 
packings remain unambiguously hyperstatic in the hard- 
sphere limit. 
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FIG. 4. Distributions of the coordination number shift to 
lower values of 2 as is increased. These results are for 
fc„ = 2 ■ lO^mg/d, e = 0.88 and <;/>' = 0.13. 
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FIG. 3. Final volume fraction 0^ versus fi for the same 
parameters as in Fig. B. 



When a static packing with a surface tilt near the an- 
gle of repose was generated by cessation of flow down an 
inclined plane similar results for z and (p were ob- 
tained. For fc„ = 2 • lO^mg/d, e = 0.88, and ^ = 0.50, 
such packings gave z = 4.69 and (p — 0.594, compared 
to z = 4.90 and 4> = 0.61 for packings presented in this 
study. 



B. The Radial Distribution Function 

The radial distribution function (RDF), g{r), for fc„ = 
2 • lif'mg/d and e = 0.88, is plotted in Fig. || for sev- 
eral values of /i. The characteristic split second peak, 
indicating short-range order out to second neighbors, is 
evident. For /i = 0, g{r) is essentially identical to that ob- 
tained for random close packing (RCP), at volume frac- 
tion 4>f « 0.64 [|l5|. As ^ increases, the secondary peaks 
in g{r) diminish, as seen in the inset. 
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FIG. 5. The radial distribution function g{r) for spheres 
with and without friction for fc„ = 2 • lO^mg/d and e = 0.88. 

The first peak of g(r) is of particular interest, since 
near-contacts with r/d just over 1 play an important role 
in the dependence of z on the stiffness of the spheres (See 
Sec. [II C). Figure |6| reveals a square-root singularity of 
the RDF near r/d — 1, i.e., 



1 



< 



«1, 



(6) 



.(i / ' d 
with a = 0.52±0.03. This singularity has apparently not 
been reported elsewhere, although we have also verified 
its presence in the RDF of hard-sphere packings provided 
to us jl^. Note that this singularity is integrable, and is 
distinct from particles actually in contact. 




FIG. 6. A logarithmic plot of g{r) vs. (r — d)/d (for 
k„ — 2 ■ Wmg/d and e = 0.88) reveals a power-law singular- 
ity with exponent a ~ 0.5 for both frictionless and frictional 
spheres. Straight line has a slope of -1/2. The same power-law 
is observed up to the largest values of fc„ studied. 



C. The Hard-Sphere Limit 

We next investigate the effect of the finite stiffness of 
the spheres on the packings, and on the average coor- 



dination number z in particular. Let us assume that 
the packings formed by these stiff elastic spheres are sta- 
tistically equivalent to packings that would be obtained 
by first forming a truly hard-sphere packing and sub- 
sequently allowing elastic relaxation. Due to the slight 
compression of the spheres of finite stiffness under grav- 
ity, we expect to see an increase in coordination number 
during this elastic relaxation. Since the typical compres- 
sive strain of a sphere under the same loading condi- 
tions scales as l/fc„, we expect that a finite fraction of 
neighbors in the hard-sphere packing that were within a 
distance d[l + 0(mg/knd)] of each other form new con- 
tacts upon elastic relaxation. The number of such near- 
contacts in the hard-sphere packing can be computed by 
integrating the hard-sphere RDF goo{'f) over the range 
r/d e (1, 1 + mg/knd), yielding an effective coordination 
number 



(^n) 



knd\ " ^ Kd^^^ 
mg J rag 



where Zao is the coordination number in the hard-sphere 
limit, a„ is a constant, and the exponent [cf. Eq.(|6)] 



a„ = 1 — a. 



(8) 



Thus, there is a power-law correction to the apparent 
coordination number, with an exponent that depends on 
the nature of the singularity in the first peak of g{r)- A 
numerical fit of the data to Eq. (^ shown in Fig. |^(a) 
results in a„ = 0.498 ± 0.002 and Zoo = 6.01 ± 0.02, in 
excellent agreement with the isostaticity hypothesis, as 
well as the exponent relation Eq. (||). Makse et al. also 
found that z approaches 6 as the stress goes to zero in 
their numerical studies of compressed spheres . 

Armed with this insight, we apply a similar analysis to 
frictional packings, with results presented in Fig. 0(b). 
Although the RDF of frictional packings appears to have 
the same square-root divergence near r = d (see Fig. |^) , 
the numerical fit to Eq. (|7|) in the presence of friction 
yields a different exponent a/ « ~l/4, resulting in a 
slower approach to the hard-sphere limit. Thus, the 
exponent identity Eq. (|8|) does not hold for frictional 
spheres. Moreover, in contrast to the frictionless case, 
the hard-sphere limits remain firmly above the isostatic 
value of four, and vary as a function of /i and e. 

Even though the rather unlikely scenario of a further 
crossover to isostaticity at extreme stiffnesses cannot be 
entirely ruled out, it must be pointed out that the stiffest 
spheres with fc„ = 2 • lO^mg/d in these packings expe- 
rience strains S/d 10~®. This should be compared to 
the strain of a "typical" grain, i.e. a glass sphere with 
a 100 micron diameter, under just its own weight on the 
Earth's surface: S/d « {pgd/ E^/^ . With the Young's 
Modulus EthQ- lO^" Pa, p w 2 • 10^ kg/m^ and g w 10 
m/s^, this strain is about 10"''. Thus, even if isostaticity 
is ultimately restored, the relevance of the hard-sphere 
limit for real granular systems is still questionable. 
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FIG. 7. Bulk averaged coordination number z as a func- 
tion of particle hardness fc„, for 0' = 0.13: (a) For frictionless 
spheres, where the extrapolation to hard-spheres implies iso- 
static packing (z=6); (b) For [/x, e] = [0.50, 0.88] (solid circles), 
[0.50, 0.50] (triangles), [0.50, 0.26] (solid squares), and [10.0, 
0.26] (solid diamonds), where the hard-sphere limit leaves the 
packings hyperstatic, with coordination numbers that depend 
on fj, and e. 



D. Plasticity of Contacts 

One potential explanation for the hyperstaticity of 
these frictional packings is the loss of degrees of free- 
dom for the tangential forces in contacts that have be- 
come "plastic" such that Ft — /iF„. If a finite fraction 
of the contacts satisfied this condition, the isostaticity 
condition would need to be modified However, as 
demonstrated by the distribution of the plasticity index 
( = of the contacts in Fig. |[ almost all contacts 
in the static packings are below their frictional threshold 
C = 1, eliminating this possibility. A similar distribution 
of C was observed for a static packing created from flow 
arrest g. 

For /X > 1, the distribution of the contact force ratio 
Ft/ Fn indeed becomes independent of /i, which manifests 
itself as a collapse of -P(C) when plotted against Ft/F„ 
(not shown). This result is in accord with the observa- 



tion in Sec. [II A that packings in this range of /i behave 
effectively the same as for systems with fi — oo. 
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FIG. 8. Probability distribution -P(C), normalized by its 
maximum value, for various values of /i. Curves for ^j. > I 
would collapse if plotted against Ft/Fn instead of (. 



E. Effect of dissipation and initial conditions 

The dependence of the packing geometry on coeffi- 
cient of restitution e is also interesting. Unlike fc„ and 
/i, changing e would not change the configuration of a 
static packing after it has stopped - it only affects the 
relaxation dynamics by increasing the removal rate of ki- 
netic energy. In this sense, changing e is like changing 
the quench rate of a supercooled liquid as it undergoes a 
glass transition. For very large quench rates, the system 
might be expected to stop immediately upon forming the 
minimum number of contacts necessary to achieve static 
mechanical equilibrium. 

Similarly to the effect of e, we find that the initial start- 
ing densities also affect the final packing. In Fig. ^ we 
plot the variation in the final packing fraction (f>f , and co- 
ordination number z, as a function of the starting density 
0' . We find that more dilute starting states lead to more 
compact final states. This behaviour may be due to the 
increase in potential energy the system receives when it 
is more dilute, converting into kinetic energy of the par- 
ticles during settling, and enabling them to explore more 
of the phase space on their way to achieving a preferred 
packing. An empirical fit to the packing fraction is given 
by, 



a/ 



0.5778 -f 0.0567 exp(-4.30' 



(9) 



which is similar to the empirical fit for model 2D systems 
proposed in Ref. [0. However, it should be noted that 
for frictional spheres, extrapolation of the coordination 
number to the limit cjf = (f)f does not result in isostaticity 
as observed in Ref. [p7[. 
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FIG. 9. Dependence of the final packing fraction (p-^ and 
coordination number z on the initial packing fraction of the 
falling particles (^', for fc„ = 2 ■ IQpmg/d, e = 0.88, and 
M = 0.50. 



In light of the dependence of the final state on such 
parameters as e and the breakdown of the expo- 
nent identity Eq.(||) for frictional packings is perhaps 
not surprising. Packings obtained by the sedimentation 
of hard-spheres followed by elastic relaxation are prob- 
ably not statistically equivalent to packings of particles 
that are elastic from the start, due to the strong history- 
dependence of the final states obtained. 

In general, for hyperstatic packings the force network 
is not uniquely determined by the packing and the load- 
ings on the particles. It follows that the determination 
of the force network for hyperstatic packings of perfectly 
rigid particles is an ill-posed problem . Thus the order 
of the limit fc„ ^ oo and the preparation of the packing 
cannot be commuted for frictional spheres. More bluntly, 
perfectly rigid sphere systems with friction are unlikely 
to capture the mechanical properties of packings of fric- 
tional, elastic spheres, even in the limit of extremely large 
rigidity of these latter particles. 



IV. CONCLUSIONS 

We have studied large scale packings of spherical grains 
of varying hardness, friction coefficient, and coefficient 
of restitution, formed by sedimentation. We accounted 
for the systematic variation with particle stiffness and 
were able to infer properties of hard-sphere packings. Al- 
though frictionless hard-spheres appear to form isostatic 
packings regardless of construction history and restitu- 
tion coefficient, frictional packings achieve a multitude of 
hyperstatic packings that depend on system parameters 
and construction history p9| . The coordination number 
reduces smoothly from z = 6 as the friction coefficient 
is increased, contrary to the hypothesis of isostaticity in 
such packings. 
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